# -*- coding: utf-8 -*-
"""
Created on Sat Oct 16 09:31:47 2021

@author: zhuo木鸟


"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate, interpolate

# 支持中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

def draw_herbs(wave_number, absorbance, no):
    '''
    根据波速和吸光度，画出某药材的中红外谱图像
    '''
    # 画出编号 no 的药材的中红外谱图像
    plt.plot(wave_number, absorbance, label=f'药材#{no}')
    #标出X轴和Y轴的标准
    plt.xlabel('波数', fontfamily='simhei', fontsize=16)
    plt.ylabel('吸光度', fontfamily='simhei', fontsize=16)
    # 画出图例
    plt.legend()

def show_herbs(wave_number, herbs_data):
    '''
    展示编号为0 100 200 400的药材的中红外谱图
    '''
    fig = plt.figure()
    for i in [0, 100, 200, 300, 400]:
        draw_herbs(wave_number, herbs_data.iloc[i, :], no=i)
    # 画图网格
    plt.grid()

    path = r'../pictures/附件1草药中红外谱展示图.png'
    # 保存图片
    plt.savefig(path)

    plt.show()

def return_distance(wave_number, y1, y2, interpolate_flag=True, step=None):
    '''
    根据积分的方法，求中红外谱曲线的距离
    
    其实也可以不用插值法来提高精度，甚至可以降低精度。。
    '''
    if step < 1:
        # 我要降低精度，提高运行速度，可以将step设为大于1，从而进行等差采用
        if interpolate_flag:
            # 为了减少运行时间，
            # 通过插值法填充数据集
            interpolate_fun_y1 = interpolate.interp1d(wave_number, y1, kind='cubic')
            interpolate_fun_y2 = interpolate.interp1d(wave_number, y2, kind='cubic')
            # 产生用于插值法的X值
            wave_number = np.arange(wave_number[0], wave_number[-1], step=step)
            y1 = interpolate_fun_y1(x)
            y2 = interpolate_fun_y2(x)
        else:
            pass
    else:
        wave_number = wave_number[::step]
        y1 = y1[::step]
        y2 = y2[::step]
    # 根据 梯形法 求解数值积分
    delta_y = np.abs(y1-y2)
    distance = integrate.trapz(delta_y, wave_number, dx=None)
        
    return distance

def initial_centroids(wave_number, datasets, k=3):
    '''
    用于初始化聚类中心
    '''
    # 防止数据在下述处理中被修改
    datasets = datasets.copy()
    centroids = np.zeros((k, datasets.shape[1]))
    # 随机抽取一个样本作为第一个初始聚类中心
    sample = datasets.sample(n=1)
    datasets.drop(index=sample.index, inplace=True)
    centroids[0, :] = sample.values
    

    for i in range(1, k):
        # 初始化权重，用以决定随机抽取样本作为聚类中心的权重
        distance_weight = np.zeros(datasets.shape[0])
        for idx in range(datasets.shape[0]): 
            # 初始化距离
            dist = 0
            for j in range(i+1):
                dist += return_distance(wave_number, datasets.iloc[idx, :], \
                                       centroids[j,:], \
                                       interpolate_flag=False, \
                                           step=3)
            distance_weight[idx] = dist**2
        
        sample = datasets.sample(n=1, weights=distance_weight)
        datasets.drop(index=sample.index, inplace=True)
        centroids[i, :] = sample.values
        
    return centroids
    
    
def Kmeans(wave_number, datasets, k=3, interpolate_flag=False, step=3):
    # 曲线的聚类
    # 初始化每个数据的所属聚类簇
    clusters = np.zeros(datasets.shape[0])
    # 采用 kmeans++ 算法，初始化聚类中心
    centroids = initial_centroids(wave_number, datasets, k=k)
    # 初始化聚类簇变化距离
    centroids_diff_flag = 1

    while centroids_diff_flag:
        min_dist = np.inf
        # 若聚类中心变化，继续更新聚类簇
        for idx in range(datasets.shape[0]):
            y = datasets.iloc[idx, :]
            for centroid_idx, centroid in enumerate(centroids):
                # 根据积分公式求取当前数据与聚类中心的距离
                dist = return_distance(wave_number, datasets.iloc[idx, :], \
                                       centroid, \
                                       interpolate_flag=False, \
                                           step=3)
                if min_dist > dist:
                    # 真心当前数据所属的聚类簇
                    min_dist = dist
                    clusters[idx] = centroid_idx
        
        # 求出均值作为新的聚类中心
        if len(set(clusters)) != k:
            # 若有一个聚类中心是多余的, 则采用上述的初始化方法，重新初始化聚类中心
            new_centroids = initial_centroids(wave_number, datasets, k=k)
        else:
            # 求个子数据的均值，作为新的聚类中心
            new_centroids = pd.DataFrame(datasets).groupby(by=clusters).mean().values
            
        # 
        if np.count_nonzero(centroids - new_centroids) == 0:
            centroids_diff_flag = 0
        else:
            centroids = new_centroids
            
    return centroids, clusters

def calculate_dist_sum(wave_number, datasets, centroids, clusters):
    '''
    求出当前，聚类簇数下数据以及各自聚类中心的距离之和
    '''
    dist_sum = 0
    for idx in range(datasets.shape[0]):
        dist = return_distance(wave_number, datasets.iloc[idx, :], \
                               centroids[int(clusters[idx])], \
                               interpolate_flag=False, \
                               step=3)
        dist_sum += dist
    return dist_sum

def elbow_method(wave_number, datasets, k_list=range(3,6)):
    '''
    画出不同的K下,数据与其各自的聚类中心的距离，从而根据拐点原则选出最佳的聚类簇数 K
    '''
    dist_sum_list = [] 
    for k in k_list:
        centroids, clusters = Kmeans(wave_number, datasets, k)
        dist_sum = calculate_dist_sum(wave_number, datasets, centroids, clusters)
        dist_sum_list.append(dist_sum)

    # 画出 elbow 图
    plt.plot(k_list, dist_sum_list)
    #标出X轴和Y轴的标准
    plt.xlabel('聚类簇数K', fontsize=16)
    plt.ylabel('数据与各自聚类中心的距离之和', fontsize=16)
    plt.title('Kmeans聚类的 Elbow 图', fontsize=16)
    plt.grid()
    plt.xticks(k_list)
    
    path = r'../pictures/Kmeans 聚类不同 k 下的 elbow 图.png'
    # 保存图片
    plt.savefig(path)
    plt.show()


def draw_herbs(wave_number, absorbance, no):
    '''
    根据波速和吸光度，画出某药材的中红外谱图像
    '''
    # 画出编号 no 的药材的中红外谱图像
    plt.plot(wave_number, absorbance, label=f'药材#{no}')
    #标出X轴和Y轴的标准
    plt.xlabel('波数', fontfamily='simhei', fontsize=16)
    plt.ylabel('吸光度', fontfamily='simhei', fontsize=16)
    # 画出图例
    plt.legend()

def draw_centroids(wave_number, centroids, clusters):
    '''
    画出聚类中心的中红外谱图
    '''
    for i in set(clusters):
        # 画出聚类中心 i 的中红外谱图像
        type_cluster = int(i)
        plt.plot(wave_number, centroids[type_cluster, :], label=f'聚类中心#{type_cluster}')
        #标出X轴和Y轴的标准
    plt.xlabel('波数', fontfamily='simhei', fontsize=16)
    plt.ylabel('吸光度', fontfamily='simhei', fontsize=16)
    plt.grid()
    # 画出图例
    plt.legend()
    path = r'../pictures/聚类中心图' + f'k={k} 时.png'
    plt.savefig(path)
    plt.show()


if __name__ == '__main__':
    # 读取附件 1 
    path = r'../datasets/附件1.xlsx'
    herbs_data = pd.read_excel(path, index_col=0)
    # 提取数据中的中红外谱的波数
    wave_number = np.array(herbs_data.columns)
    
    # 展示编号为0 100 200 300 400的药材的中红外谱图
    show_herbs(wave_number, herbs_data)
    
    # 画 elbow 图，用于找出最佳的聚类簇数K
    elbow_method(wave_number, herbs_data, k_list=range(3,21))
    
    
    # 修改这里的K，可以求解不同聚类簇塑像的情况
    k = 8
    centroids, clusters = Kmeans(wave_number, herbs_data, k=k)
    
    # 画出聚类中心
    draw_centroids(wave_number, centroids, clusters)
    
    herbs_type = pd.DataFrame(data=clusters, index=herbs_data.index)
    path = r'../results/附件1各药材编号对应的药材种类' + f'k={k} 时.xlsx'
    herbs_type.to_excel(path,  header=['药材种类'])
